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The kinetic behavior of a phase field model of electrochemistry is explored for advancing (elec- 
trodeposition) and receding (electrodissolution) conditions in one dimension. We described the 
equilibrium behavior of this model in [J. E. Guyer, W. J. Boettinger, J. A. Warren, and G. B. Mc- 
Fadden, "Phase field modeling of electrochemistry I: Equilibrium", cond-mat/0308173]. We examine 
the relationship between the parameters of the phase field method and the more typical parameters 
of electrochemistry. We demonstrate ohmic conduction in the electrode and ionic conduction in the 
electrolyte. We find that, despite making simple, linear dynamic postulates, we obtain the nonlinear 
relationship between current and overpotential predicted by the classical "Butler- Volmer" equation 
and observed in electrochemical experiments. The charge distribution in the interfacial double layer 
changes with the passage of current and, at sufficiently high currents, we find that the diffusion 
limited deposition of a more noble cation leads to alloy deposition with less noble species. 

PACS numbers: 81.15.Aa, 81.15.Pq, 82.20.Wt, 82.45.Qr 



I. INTRODUCTION 

In Ref. [1], we developed an equilibrium phase field 
model of an electrochemical system. In this paper, we 
examine the dynamic aspects of that model. Models of 
phase transformations can be broadly categorized into 
sharp or diffuse interface approaches. Sharp interface 
models treat the transition between phases as mathe- 
matically abrupt. Diffuse interface models assume that 
the phase interface has some finite thickness over which 
material properties vary smoothly. Both cases are sim- 
plifications of the physical interface between phases, in 
which properties vary over some finite, atomic-scale dis- 
tance which is often smaller than assumed in diffuse in- 
terface models. Traditional equilibrium models of elec- 
trochemical interfaces take the interface between phases 
(the transition between "electrode" and "electrolyte" ) to 
be abrupt, but frequently consider the distribution of 
charge and electrostatic potential to be diffuse in the 
electrolyte, as by the Gouy-Chapman-Stern model [2]. 
Dynamic models of electrochemistry typically ignore the 
details of the charge distribution at the interface and em- 
ploy a fully sharp model that we summarize in Section II. 

The phase field technique is one particular diffuse in- 
terface approach. The method employs a phase field vari- 
able, which is a function of position and time, to describe 
whether the material is one phase or another, e.g., solid 
or liquid. The behavior of this variable is governed by 
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a partial differential equation (PDE) that is coupled to 
the relevant transport equations for the material. The 
interface between the phases is described by smooth but 
highly localized changes of this variable. This approach 
avoids the mathematically difficult problem of applying 
boundary conditions at an interface whose location is 
part of the unknown solution. The phase field method 
is powerful because it can easily treat complex interface 
shapes and topology changes. Our long term goal is to 
treat the complex geometry, including void formation, 
that occurs during electrodeposition in vias and trenches 
for on-chip metallization and in the dendritic structures 
that form during battery recharging. Phase field meth- 
ods will allow the rigorous examination of the interplay 
between current, potential gradients, curvature, and ad- 
sorption in intricate geometries. 

There is a rich body of literature of sharp interface 
models of electrodeposition, which we will sketch briefly 
in Section II, but the application of diffuse interface tech- 
niques to the motion of electrochemical interfaces has 
been relatively limited. Dussault and Powell have applied 
phase field techniques to the modeling of electrochemical 
processes in steel slags [3, 4], but their approach neglects 
the effects of charge at the interfacial double layer. As a 
result, they are able to model much larger domains and 
much longer time spans than we present here, but the 
essential physics of the electrocapillary interface are not 
examined. Wheeler, Josell, and Moffat have performed 
a level set analysis of so-called "superconformal" elec- 
trodeposition in high aspect ratio features, with particu- 
lar emphasis on the role of additives [5] . Like phase field 
models, level set techniques allow the treatment of com- 
plex morphologies, such as the formation of voids dur- 
ing trench filling, but the motion of the interface is han- 
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died phcnomcnologically rather than physically, as by the 
phase field approach, so again the structure of the dou- 
ble layer is not considered. Bernard, Plapp, and Gouyet 
have recently presented a lattice-gas model of an electro- 
chemical system [6, 7] that exhibits many of the same 
intcrfacial and dynamic behaviors that we find in this 
paper, as well as exhibiting the early stages of dendritic 
growth. This type of discrete modeling should provide a 
useful bridge between an atomistic view of the electro- 
chemical interface and the continuum approach of phase 
field models. 

To place our results in context, Section II outlines the 
traditional sharp interface description of electrodeposi- 
tion. Section III presents the dynamic postulates govern- 
ing the evolution of the phase, concentration, and elec- 
trostatic potential fields which we proposed in Ref. [1]. 
Section IV describes our numerical approach and bound- 
ary conditions. Section V discusses the selection of mate- 
rials parameters, including the relationship between the 
phase field mobility and the Butler- Volmer exchange cur- 
rent of traditional electrochemical modeling. Section VI 
presents the results of numerical calculations in one spa- 
tial dimension that span a range of electrodeposition and 
electrodissolution conditions. 



II. SHARP INTERFACE APPROACH FOR 
ELECTRODEPOSITION 

In Ref. [1], we present a phase field model of the equi- 
librium between an electrode phase a and an electrolyte 
phase P, consisting of a set of four charged components, 
c", M +m , N +n , and A~ a . A superscript a denotes that 
the quantity is evaluated in the bulk electrode (metal) 
phase and a superscript (3 denotes that the quantity is 
evaluated in the bulk electrolyte phase. At equilibrium, 
the difference in potential between the electrode and 
the electrolyte in an n-component system is given by 
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where A[i° — fi° a — \x? and ji° is the chemical potential 
of pure component j in the respective phase, Zj is the 
valence of component j, Xj is the mole fraction of com- 
ponent j, T is Faraday's constant, R is the molar gas 
constant, and T is temperature. Eq. (1) is the general- 
ization, for all of the components, of the Nernst equation 
of traditional electrochemical analysis. The equation is 
normally only written for the electroactive species. In 
Ref. [1] we explain the origin of an equation for each 
component in the system and the relationship between 
the term proportional to A[i° and the standard cell po- 
tential. 

When current is passed through the interface, the po- 
tential difference Acj) shifts. Alternatively, when a po- 
tential difference other than the equilibrium value is im- 
posed, current will pass and the interface will move. The 



shift in the potential difference across the interface (ex- 
cluding the ohmic drops across the bulk phases) is re- 
ferred to as the overpotential 77 [8, 9]. 

For an electrochemical system with only one mono- 
valent electroactive species M + , chemical reaction rate 
theory gives the relationship between current density i 
and total overpotential rj for a planar interface as [8] 
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The first term in the curly brackets represents the an- 
odic/oxidizing reaction and the second term represents 
the cathodic/reducing reaction. C M + is the concen- 
tration of cations M + in the electrolyte at electrode- 
electrolyte interface and is the bulk electrolyte con- 
centration of cations M + at the edge of the diffusion 
boundary layer. The exchange current density iq and the 
transfer coefficient v characterize the facility and sym- 
metry of the forward and reverse reactions. The current 
density i = i • n, where the normal n points from a into f3 
and i is the current density vector. Thus, positive values 
of i result in dissolution. 

If the diffusion field can be assumed to be linear, the 
implicit dependence of C M + on i can be eliminated in 
Eq. (2), giving 
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This expression can be rearranged to give i as an explicit 
function of 77, which is useful for comparison to our phase 
field results. Linearity of the concentration profile is ap- 
propriate only if the interface velocity is much less than 
D m +/5d, where D M + is the diffusivity of M + and 5d is 
the thickness of the diffusion boundary layer. The limit- 
ing deposition current in m is determined by the complete 
depletion of M + in the electrolyte at the interface, such 
that 
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The classical "Butler- Volmer" equation of electrochem- 
istry is a special case of Eq. (3) in which the effects of 
mass transfer are neglected (i/ium —> 0). 

For small ovcrpotcntials, the linearized form of Eq. (3) 

is 
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We will use this relationship in Section V C to relate 
io to the parameters of our phase field model. When 
\rj\T/RT > 1 and i <C ium, Eq. (4) reduces to 
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The quantities (1 — v)T jRT and —vTjRT are known as 
the anodic and cathodic "Tafel slopes" from the slopes 
of the lines when ln|i| is plotted against r\. These slopes 
can be used to deduce experimental values for v. 

Eq. (2) was originally derived from reaction rate theory 
to explain experimentally observed current-overpotential 
behavior. More recently, atomistic and quantum me- 
chanical treatments of electron and ion transfer reactions 
have been performed to replace this chemical reaction 
rate approach [10]. These treatments have led to a bet- 
ter physical understanding of the phenomenological con- 
stants v and «0j but they do not fundamentally alter the 
form of Eqs. (2) and (3). 
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Poisson's equation 
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must also be satisfied everywhere, where the charge den- 
sity is 



P = F^ZjCj. 



(12) 



The mobilities Mj and M$ will be related to the param- 
eters of electrokinetics in Sections VB and VC. 



III. MODEL 
A. General Kinetic Equations 

In Rcf. [I], we performed a variational analysis to de- 
rive the governing equations for the equilibrium electro- 
chemical interface. We also postulated the simplest time 
dependent forms of those governing equations that guar- 
antee a decrease in total free energy with time t. We 
restate those dynamic postulates here. The time varia- 
tion of the phase field £ is given by 
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where fy is the Helmholtz free energy density per unit 
volume, k,£ is the phase field gradient energy coefficient, 
e(£) is the dielectric constant, which we take to depend 
explicitly on the phase; because all of the fields are cou- 
pled, it will also depend implicitly on the electrolyte con- 
centration, is the mobility of the phase field. Under 
the assumption that all nonzero partial molar volumes 
are identical, the flux Jj of each component j is 
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where Jij is the electrochemical potential of species j and 
Vj is the partial molar volume of species j. We divide 
the components into electrons e~ with j = 1, which have 
V e - = 0, and substitutional species with j > 1, which 
all have the same Vj = V s =0. One consequence of this 
assumption is that X)J=2 = = constant, where 
Cj is the concentration of species j. A specific choice 
is made of a substitutional component n with nonzero 
partial molar volume to be called the reference species. 
The quantity Mj is the mobility of component j. Since 
conservation of species requires 



at 



-v- J,- 



j 



. n 



B. Form of the Dynamic Equations for Ideal 
Solution Behavior 



For simplicity, we assumed in Ref. [1] that the chemi- 
cal part of the Helmholtz free energy per unit volume is 
described by an interpolation between two ideal solutions 
of the components, 
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where the molar volume V m = (Y^=i^j) We use 
an interpolating function p (£) = £ 3 (6£ 2 — 15£ + 10) to 
bridge between the descriptions of the two bulk phases 
and a double-well function g (£) = £ 2 (1 — £) 2 with a bar- 
rier height Wj for each component j to establish the 
metal/electrolyte interface [11]. The polynomials are 
chosen to have the properties that p(Q) = 0, p(l) = 1, 
p'(0) = p'(l) = 0, and g'(0) = g'(l) = 0. The clas- 
sical chemical potential is given by fij = dfy/dCj and 
the corresponding classical electrochemical potential is 
fij = ft j + ZjT4>. 

Substituting Eq. (13) into Eqs. (7) and (8), we obtain 
the governing equation for evolution of the phase field 
under ideal solution thermodynamics 
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and the flux in the diffusion equation Eq. (9) is given by 
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TABLE I: Numeric values of the potential-independent por- 
tion of the chemical potential differences A/i°. 



hx{X?°/X?°) 


e 


-20.03 


M+ 


-3.912 


N+ 


20.01 


A" 


20.03 



MATERIAL PARAMETERS 



A. Equilibrium Material Parameters 



The flux of substitutional species does not explicitly de- 
pend on the electron concentration and the flux of elec- 
trons does not explicitly depend on the concentration 
of substitional species; the flux of substitutional species 
is affected by the displacement of other substitutional 
species, but electrons can move without displacing other 
ions. The fluxes of all species are coupled indirectly 
through the total charge distribution and Eq. (11). 



IV. NUMERICAL METHODS 

The 1-D form of the governing equations was trans- 
formed to a frame moving at a velocity v. Simulations 
were performed in a domain of length L with an initially 
abrupt interface between the bulk electrode and elec- 
trolyte phases at x = L/2, such that £,\ x< l/2 = = 1 
and £U>l/2 = = 0. After choosing an initial bulk 

value for C^ + ; the remaining initial bulk Cf and Cf 
were the equilibrium values obtained by equating the 
bulk electrochemical potentials p,j [1]. The boundary 
condition on the phase field is n • V£ = at both ends 
of the solution domain. At the electrolyte end, we set 
cf> = and at the electrode end we specify i. At the lead- 
ing edge of the moving frame, we model the stirred bulk 
electrolyte by applying a fixed concentration boundary 
condition. At the trailing edge of the frame, we discard 
the material leaving the frame by setting the divergence 
of the species fluxes to zero. 

Equations (9), (11), (14), and (15) were solved with 
explicit finite differences. Spatial derivatives were taken 
to second order on a uniform mesh. Transient solutions 
were integrated numerically with an adaptive, fifth-order 
Rungc-Kutta time stepper (based on odeint of Ref. [12]) 
until a steady state was achieved (current became con- 
stant). We have defined steady state in our simulations 
as the point when each Jj — vCj were uniform to within 
0.1%. Because v is an unknown result of the simulation, 
the frame velocity was adjusted at each iteration to keep 
the interface stationary in the frame. 



We examine the dynamic behavior of a four component 
model under a different set of thermodynamic parameters 
than described in Ref. [1]. In this paper, all components 
have valence Zj = ±1. We are primarily interested in the 
electrodeposition of the more noble cation M + , where the 
less noble cation N + and the anion A~ make up the bulk 
of the supporting electrolyte. This electrolyte contain- 
ing only charged species represents a molten salt system. 
The presence of the second cation N + introduces the pos- 
sibility of alloy deposition. 

We take the partial molar volume of the "substi- 
tutional" components (M + , N + , and A - ) as V s = 
10" 5 m 3 /mol. Equation (1) states that for any given X" 

and Xj, there is some potential difference A<fi between 
that bulk phases that is related to the chemical poten- 
tial difference of the pure components A/i°. Conversely, 
we showed in Ref. [1] that we can establish a value for 
A/ij if we know A(j> for some particular X" and X? , for 
instance the standard state values 
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The voltage-independent portion of A/i° is given in Ta- 
ble I. In this paper, we take Ac/) to be zero. Following 
Ref. [1], this implies that the equilibrium state for this 
material system at the standard state concentration is 
near the point of zero charge. The mole fraction ratios 
in Table I of the normally-electroinactive species are cho- 
sen to give the corresponding small standard state mole 
fractions as X ° = X°° + = X"° = 10" 9 . 

c N + A 

To permit a convenient graphical display of bulk equi- 
librium, we invoke charge neutrality to transform the four 
components {e~, M + , N + , A~} into an alternate set of 
four components that are charge- neutral, {M, N, MA, 
NA}. Wc plot the equilibrium phase diagram in terms of 
these transformed components in Figure 1. Equilibrium 
states exist only between -0.5138 V < Acj) < +0.1005 V. 
It can be seen that over the majority of the potential 
range, from -0.4 V < A<fi < V, that the equilibrium 
is between an electrode of essentially pure M and a NA 
electrolyte containing a dilute concentration of MA. At 
the positive A<fi extreme, the equilibrium is between M 
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electrolyte 



MA 




electrode 

FIG. 1: Potential-composition phase diagram for the param- 
eters in Table I, illustrating the bulk equilibrium between 
a M electrode and an electrolyte containing MA salt dis- 
solved in NA. Tie-lines denote different values of the quan- 
tity (Acj) — Acj)°). The inset shows the position of this charge 
neutral phase diagram within the quaternary domain of the 
charged species. 



TABLE II: Parameters characterizing the equilibrium inter- 
face, eo is the permittivity of free space. 



parameter 


value 




7.2 x 10" 11 J/m 


Wj£2...n 


3.6 x 10 5 J/mol 


W e - 


J/mol 


C 


8e 



and MA and at the negative A</> extreme, the equilib- 
rium is between a phase of M and N and a solution of N 
and NA for this choice of X° . 

Table II lists our choice of the parameters that char- 
acterize the thickness and energy of the electrode- 
electrolyte interface. Our assumption that the barrier 
heights Wj are equal for the substitutional species and 
zero for electrons is discussed in Ref. [1]. 



B. Single Phase Transport Properties (values for 

To identify the mobilities Mj, we examine single-phase 
systems. In a single-phase electrode, £ = p(£) = 1. In 
a single-phase electrolyte, £ = = 0. In either phase 
= Vp(£) = Vg(£) = 0. We thus can write the 



fluxes in Eq. (15) as 
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The flux of component n balances the other fluxes such 
that 
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We first consider an electrolyte with = 0. If we 
compare the resulting form of Eq. (17a) with the classical 
diffusive flux equation with diffusivitics D^, 
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IS 

D.j,jC n Cj 



elements of as 



Mj = 



RT(C n +Cj) 



J 



.n-l 



(21) 



For simplicity, we assume the diagonal elements of Dij 
are constants, thus inducing a concentration dependence 
in the mobilities as defined by Eq. (21) and in the off- 
diagonal Dij's. 

We next consider an electrode with all VCj = 0, where 
the current is entirely carried by the electromigration of 
electrons. The resulting form of Eq. (17b) 



J e - = -M e -z e -FV4> 
can be substituted into Eq. (18) to give 
i w -z^T 2 M c -V(j) 
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By comparison with Ohm's law, i = —aV<f>, we readily 
see that the electron mobility 
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Thus Eqs. (21) and (24) relate the M/s to the electronic 
conductivity and ionic diffusivities. 
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On substitution of Eq. (21) into Eq. (17a), we see that 
the electromigration flux (due to gradients in (f>) within 
the electrolyte is 
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This is just as expected from traditional electrochemical 
theory, in the dilute limit where C„/(C„ + Cj) « 1. We 
will find in Section VIA that, for our supported ionic 
electrolyte and our electronic conducting electrode, we 
are justified in neglecting the contributions of the elec- 
tromigration current in the bulk electrolyte and of the 
diffusion current in the bulk electrode. 

It is interesting to note that the conductivity predicted 
by Eq. (24) is completely analogous to that predicted by 
the Drudc model (and by the Fermi-Dirac model, for that 
matter) [13] 
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where m e - is the mass of the electron. The relaxation 
time r can only be determined by quantum mechanical 
means and is simply an unknown constant in classical 
models of electron transport. Following an analysis for 
the electrons similar to that which gave us Eq. (21), we 
find that we can describe the mobility of electrons M e - 
in terms of a constant diffusivity of electrons D e - 
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In a single-phase conductor with uniform concentrations, 
(l+t^jC e -) is a dimcnsionless constant of order 1. We can 
see that D c - / (RT) is dimensionally equivalent to r/m e - 
and all other terms in Eqs. (24) and (26) are identical. 
The room temperature conductivity of silver of approx- 
imately 6 x 10 7 0" 1 m" 1 results in D c - « 8 x 10" 5 m 2 /s 
and M e - w 6 x 10" 3 mol 2 /(J sm). We observe that one 
of the weaknesses of the Drude model is that it fails to 
predict the a ~ T" 1 dependence found in experiments 
without making some unsatisfactory ad hoc assumptions; 
this dependence arises naturally in our fundamentally 
thermodynamic formulation. 



C. Interfacial Kinetics (value for M^) 

Along with the transfer coefficient v, the exchange cur- 
rent io characterizes the kinetics of the interface and we 
hypothesize that it has an intimate relationship to the 
phase field mobility Mg. To test this hypothesis for our 
model, we examine Eq. (5) and plot rj obtained from 
steady-state calculations against M^ 1 for various Djj 
and small values of i in Figure 2. If we scale length by 



L (the length of the solution domain), time by L 2 /Djj, 
energy density by RT/V S , and potential by RT j T , we 
find that all of the points satisfy the linear relationship 

t/*2 RTV L 

rj/i = (6.610 ±0.QQQ)Mf 1 -^-- + (4563 ± 1) 

Lij~ 

Comparison with Eq. (5) reveals that 
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and 



iiim = (2.191 ± 0.0004) x 10" 4 ^i^. 

LV S 

w 2.11 x 10 6 A/m 2 when D j3 = 10" 9 m 2 /s 

(30) 

Eq. (29) confirms our hypothesis that io is directly re- 
lated to M^. Comparing Eq. (30) to Eq. (4), and taking 
= 10 mol/m 3 , we see that this implies that the diffu- 
sion boundary layer thickness is 8d = (0.4564±0.0001)L. 
This is very close to the thickness of the electrolyte, which 
validates that we are computing the diffusion field cor- 
rectly (because we are modeling a diffuse interface, the 
electrolyte thickness is somewhat less than 0.5L). The 
thinness of the diffusion boundary layer in our calcula- 
tions gives rise to a limiting current that is much larger 
than encountered in physical systems, but the mechanism 
is the same. 

Table III displays the kinetic parameters of the phase 
field model and typical values of the corresponding physi- 
cal quantities. If physical values are used for some kinetic 
parameters, then the computation time is too long, so the 
values used for our numeric simulations are also listed. 



VI. NUMERICAL RESULTS AND DISCUSSION 

Our purpose is to show consistency of behavior with 
sharp interface models of electrochemical systems so that 
future 2-D and 3-D computations treating more complex 
phenomena can be performed with confidence. In this 
section, wc examine the behavior of our model in the 
bulk phases, explore the current-overpotential behavior, 
and demonstrate the electrodeposition of alloys at high 
applied currents. 

The interfacial region of a representative steady-state 
solution, with i = -100 A/m 2 , is displayed in Figure 3. 
The phase field £, concentrations Cj, charge density p, 
and electrostatic potential 4> are plotted against the same 
x-axis. The velocity of the moving frame is indicated 
with a marker on the £ curve at £ = 0.5. To highlight 
the location of the interface, <?(£) is mapped onto the 
background in gray. We can see that the concentrations 
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FIG. 2: (a) Dimensional and (b) dimensionless relationship between overpotential r\ and current i as a function of the inverse 
phase field mobility M% 1 when = 10 mol/m 3 and Djj = 10" 9 m 2 /s. Points are plotted for each permutation of i = 

(-500, -100, 100) A/m 2 and M £ = (10~ 3 , 1.5 x 10" 3 , 3 x 10~ 3 , 10" 2 , 10"\ 1) m 3 /(Js). Points are also plotted for i = 100 A/m 2 
and M ( = (1.5 x 10" 3 , 3 X 10~ 3 , 10~ 2 , 1) m 3 /(J s) with Djj = 10" 10 m 2 /s. The line in Figure (b) is a fit to r//i = aM^ 1 + b 
with a = (6.610 ± 0.006)(V s 2 /i^ 2 ) and b = (4564 ± 1 ) (RTV S L/Djj T 2 ) . The points in the dotted box contribute to Figure 5. 



TABLE III: Correspondence between kinetic parameters used in this phase field model and those measured in experiments or 
typical of sharp-interface models. Typical physical values [9] are compared with the values used in our numerical calculations. 
The diffusivities are given for the electrolyte phase; diffusivities in the solid electrode are expected to be many orders of 
magnitude smaller. For the calculations in this paper, we treat the diagonal diffusivities as constant and uniform. To simplify 
the notation, we take Dj = Djj. No D A - is necessary because A - is the reference species in our calculations. 



phase field 


"physical" 


"numeric" 


D e - = 10" a m 2 /s 


a = 6 x 10 Y Q- L m L 


a = 750 O" 1 !!!" 1 


D M + = IO" 9 m 2 /s 


D M + = 10" 9 m 2 /s 


D M+ = 10" 9 m 2 /s 


D N+ = 10" 9 m 2 /s 


D N+ = 10" 9 m 2 /s 


D N + = 10" 9 m 2 /s 


M e = IO" 2 m 3 /(J s) 


io = (10- 16 to lO" 2 ) A/m 2 


io = 3.7 x 10 6 A/m 2 



deviate from their bulk values in a region of approxi- 
mately the same thickness as the phase field transition. 
As a result, the charged "double layer" is confined to 
this same region. The surface of the electrode has ex- 
cess e~ , whereas the surface of the electrolyte is an es- 
sentially charge-neutral NA salt with a dilute concentra- 
tion of MA. All of the species except M + are excluded 



from the region of intermediate £, giving rise to a layer 
of M + that has neither e~ nor A~ to balance the charge. 
This charge distribution gives rise to the potential step 
of approximately 0.12 V between the two phases, which 
is the expected Nernst potential of an electrolyte with 
qfL = 10 mol/m 3 . 
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FIG. 3: Interface profiles for steady state electrodeposition 
with i = -10 2 A/m 2 . The concentration profiles for N + and 
A - are almost coincident on this scale. is mapped onto 
the background in gray to indicate the location of the phase 
field interface. 



A. Fluxes 

The relative contributions of the flux due to diffu- 
sion (dependent on all the VQ) and the flux due 

to electromigration (proportional to V</>) can be dis- 
tinguished using Eq. (17). For i = -100 A/m 2 , the partial 
fluxes in the bulk electrolyte are listed in Table IV and 
those in the bulk electrode are listed in Table V. As the 
designated reference species, the flux of A~ always ad- 
justs such that the sum of the fluxes of the substitutional 
species is zero. In both phases, the concentration gradi- 
ents of M + and N + are approximately equal and opposite 
in sign to maintain charge neutrality (the concentration 
gradients of A~ and c~ are small). The diffusive fluxes of 
M + and N + are not equal and opposite in sign. The "off- 
diagonal" term for the N + flux in the electrode and for 
both the M + and N + fluxes in the electrolyte contribute 
significantly. 

Because we consider a supported electrolyte (the to- 
tal ion density is high) , V</> is small and electromigration 
does not contribute significantly to the current in the 
electrolyte. Both the magnitude and gradient of C c - are 
small in the electrolyte, such that e~ do not carry any 
significant current in the electrolyte. The current due to 
the N + flux is cancelled by that due to the A - flux, such 
that essentially all of the current in the electrolyte is car- 
ried by the diffusion of M + . In the electrode, the partial 
fluxes of the substitutional components are numerically 
zero. The concentration gradient of e _ is small in the 
electrode, giving a small diffusive flux. The bulk of the 
current in the electrode is carried by electromigration of 
e - , consistent with Ohm's law. These observations that 
the current in the electrolyte is carried by diffusion of M + 
and the current in the electrode is carried by the electro- 
migration of e _ are consistent with the approximations 
we made for the bulk phases in Section VB; i.e., bulk 
behavior is obtained at a distance of 0.5 nm from the 
interface. 



In Figure 4 we plot the profile of M + in the electrolyte, 
showing the depletion due to electrodeposition and the 
enrichment due to electrodissolution. At the highest cur- 
rent in Figure 4(a), we can see that C M + near the sur- 
face of the electrode is depleted practically to zero, giving 
rise to the limiting current behavior of Section II. The 
diffusion layer thickness o~q = 0.456L, calculated in Sec- 
tion VC, is indicated for comparison. Over the range of 
applied currents examined, the enrichment of M + during 
electrodissolution is not similarly constrained. 



C. Current-Overpotential Relationship 

In Section V C, we found that the relationship between 
current i and overpotential rj in our calculations is satis- 
fied by the linear relationship (5) when i and r\ are small. 
Now we plot i vs. r\ over a larger range of applied currents 
in Figure 5 as open squares. Equation (3) considers only 
the electroactive species, so the filled circles in Figure 5 
show the current carried by the electroactive cation i M + . 
The relationship between i M + and rj is not linear. At 
large, negative values of ry, we observe a limiting current, 
whereas for large positive values of rj, no such limiting 
current is observed and «m+ appears exponentially de- 
pendent on rj. We fit Eq. (3) to the calculated values 
of i M + and we find that i = (3.80 ± 0.08) x 10 6 A/m 2 , 
i hm = (-2.15 ± 0.06) x 10 6 A/m 2 , and v = 0.777±0.002. 
These values of iq and i\ im are within 5% of the values 
found in the linear analysis of Section VC. Because %q is 
of the same order as in m in our calculations, we do not 
observe an obvious "Tafcl slope" during electrodeposi- 
tion. Nonetheless, the transition between low current and 
diffusion-limited current cannot be fit except by the full 
form of Eq. (3) . From these results, we see that despite 
postulating a linear evolution equation for the phase field 
(Eq. (7)), we obtain the nonlinear current-overpotential 
behavior predicted by sharp-interface theories and ob- 
served in electrochemical experiments. 

The transfer coefficient v characterizes the symmetry 
of the energy barrier between the electrode and elec- 
trolyte phases. A value of v = 0.5 would mean the energy 
barrier is symmetric and that a given change in potential 
would cause the barrier to electrodeposition to change by 
the same magnitude as the barrier to electrodissolution. 
Our observed value of v = 0.78 indicates that the barrier- 
to electrodeposition is more sensitive to changes in poten- 
tial than is the barrier to electrodissolution. Although we 
do not know the functional relationship between v and 
the parameters of the phase field model, we can surmise 
that it is related to the height Wj and shape g(£) of the 
interfacial energy barriers. This will be investigated in 
the future. 

Since the exchange current is equal to the balanced 
anodic and cathodic current passed at equilibrium, it can 
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TABLE IV: Partial fluxes in the bulk electrolyte for i = -100 A/m 2 (electrodeposition) . V<£ = 6.87 x 10" 4 V/m. 



3 


C,/ (mol/m 3 ) 


VCj /(mol/m 4 ) 


Jf /(molm- 2 ^ 1 ) 


^/(molm- 2 ^ 1 ) 


Jf tal /(molm- 2 s" 1 ) 


e~ 


1.00 x 10^ 


-1.65 x 10 3 


1.65 x 10" B 





1.65 X 10 -B 


M+ 


1.00 x 10 1 


1.03 x 10 6 


-1.03 x 10 -3 





-1.03 x 10" 3 


N+ 


5.00 x 10 4 


-1.03 x 10 6 


5.18 x 10" 4 


-1.34 x 10" 6 


5.16 x 10" 4 


A" 


5.00 x 10 4 


8.28 x 10 2 






5.17 x 10" 4 



TABLE V: Partial fluxes in the bulk electrode for i = -100 A/m 2 (electrodeposition). \7cj> = 0.133 V/m. 



3 


^■/(mol/m 3 ) 


VCj /(mol/m 4 ) 


Jf /(molm-V 1 ) 


^/(molm- 2 ^ 1 ) 


Jf tal /(molm- 2 s" 1 ) 


e~ 


1.00 x 10 5 


4.97 x 10 J 


-4.97 x 10"' 


1.03 x 10" 3 


1.03 x 10" 3 


M+ 


1.00 x 10 5 


-5.93 x 10 5 


1.85 x 10" 7 





1.85 x 10" 7 


N+ 


2.04 x 10" 2 


5.93 x 10 5 


1.24 x 10" 7 





1.24 x 10" 7 


A- 


2.10 x 10" 6 


1.88 x 10 2 






-3.10 x 10" 7 



be shown that [8, 9] 



D. Alloy Electrodeposition 



(31) 



Cq is the concentration of the oxidized clcctroac- 



tive species in the bulk electrolyte, which is — 
10 mol/m 3 in our notation. C|? is the concentration of 
the reduced electroactive species in the bulk electrode, 
which is ~ l/V s in our notation. The only terms 

we cannot directly identify in our phase field model are 
the dimcnsionlcss transfer coefficient v and the rate con- 
stant ko. Noting that we found io oc in Section VC, 
from a dimensional analysis, one may expect that 



(32) 



The surface free energy found in our paper on the equi- 
librium electrochemical interface [ll is 



7 : 



dx. 



(33) 



From numerical calculations on the system in this paper 
when i = 0, we obtain a value of 7 = 0.46 J/m 2 . If we 
assume that fco in Eq. (32) is not just proportional to but 
equal to M^7 and substitute this value of the surface free 
energy and Eq. (29) into Eq. (31), we obtain v w 0.73. 
If we assume instead that the surface free energy is that 
found in models of single component solidification [14] 



7 = 



1 K^W 

18V S 



= 0.38 J/m 2 



(34) 



We examine electrodeposition of alloys by increasing 
the applied current by five orders of magnitude from 
-10 2 A/m 2 to -10 7 A/m 2 , starting from the steady state 
result of Figure 3. The fields in the vicinity of the in- 
terface are displayed at four different times in Figure 6. 
We have added a small concentration inset to each frame 
to highlight the behavior of M + in the electrolyte and a 
bar of color that represents the overall composition of the 
system. The initial potential drop across the interface of 
A<p = 0.118 V is within 2 pV of the Nernst potential for 



c 



18 

M+ 



10 mol/m 3 . At 10 ns after the step in current, 



C M + has depleted at the interface to approximately half 
its bulk value and N + has begun to accumulate at the 
electrode surface. At 200 ns, C M + has depleted essen- 
tially to zero at the electrode surface, giving rise to the 
limiting current of M + through the electrolyte. This M + 
current of approximately —2.1 x 10 6 A/m 2 is not ade- 
quate to meet the applied current of -10 7 A/m 2 . The 
surface of the electrode becomes covered with a layer very 
rich in N + and an alloy of M and N begins to deposit on 
the electrode. By 750 ns, the interfacial structure estab- 
lished at 200 ns is essentially unchanged and the original, 
pure M electrode has been completely swept from view, 
replaced by a MN alloy. 

In Figure 7, we plot the steady-state concentration of 
N + in the electrode as a function of r\. For small over- 
potentials, up to 77 ~ -0.17 V, the electrode is essentially 
pure M. At large magnitudes of 77, the fraction of N 
grows in an apparently linear fashion. 



we find that v w 0.75. In either case, this value of v is 
very close to that obtained by comparing our results to 
the the sharp interface equation (3), and is not strongly 
sensitive to the choice of 7. Although v is usually as- 
sumed to be 1/2 when no other information is available, 
it can take on any value between and 1 for an ion trans- 
fer reaction [10]. 



E. Interface Structure 

The concentration and charge distributions at the in- 
terface are sensitive to the electrodeposition conditions 
at all overpotentials or applied currents, but can be seen 
clearly in Figure 6. At % = -10 2 A/m 2 , C M +, C N + and 
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FIG. 4: Concentration of M + as a function of position in the electrolyte for different total (a) electrodeposition and (b) 
electrodissolution currents. The concentration at the interface exhibits a Nernstian shift in concentration with overpotential as 
the current is changed. The diffusion boundary layer is clearly linear over the small simulation domain. g(£) is mapped onto 
the background in gray to indicate the location of the phase field interface. The dashed vertical line indicates the thickness of 
the diffusion boundary layer S D calculated in Section VC. The concentration gradient at i — -10 7 A/m 2 gives rise to the M + 
limiting current of iii m ~ —2 x 10 6 A/m 2 ; the majority of the current is carried by other species. 



C A - in the electrolyte remain very close to their bulk val- 
ues, all the way into the interfacial region. The charge 
distribution consists of a dipolc on the electrode side, 
with very small net negative charge, and a correspond- 
ing positive charge on the electrolyte. At i = -10 7 A/m 2 , 
C M + is depleted nearly to zero at the interface and N + 
displaces essentially all of the A~ at the interface. The 
density of e~ at the surface of the electrode is much larger 
than at the lower current and the charge distribution has 
shifted to a predominantly negative charge on the elec- 
trode and a positive charge on the electrolyte. These 
changes in the charge distribution are directly tied to 
the change in overpotential, through Eq. (11). 



VII. CONCLUSIONS 

Previously [1], we developed a phase field model of the 
electrochemical interface. We performed numerical cal- 
culations on a model system like an aqueous electrolyte, 
in which the majority species in the electrolyte had no 
charge. We demonstrated that, even with a simple ideal 
solution thermodynamic description, our model exhib- 
ited charged double layer behavior, an "electrocapillary" 
relationship between surface free energy and electrostatic 



potential difference across the interface, and differential 
capacitance curves that are strongly reminiscent of ex- 
perimental measurements. 

In this paper, we have applied the same phase field 
model to electrodeposition and electrodissolution condi- 
tions. We have performed numerical calculations on a 
model system like a molten salt, with four species which 
all carry charge. We have shown: 

• the relationship between the parameters of the 
phase field model and the physical parameters of 
an electrochemical system, 

• that our model electrode carries current by elec- 
tromigration of electrons and that our model elec- 
trolyte carries current by diffusion of cations, 

• that the diffusion field in the electrolyte is essen- 
tially linear and that limiting current behavior re- 
sults, 

• that despite making linear postulates for the 
time-dependent governing equations, the current - 
overpotential relationship is non-linear and agrees 
very well with the classic sharp-interface relation- 
ship ( "Butler- Volmer" with mass transport effects) , 
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FIG. 5: Magnitude of the current i as a function of the overpotential r\, plotted on (a) linear and (b) log-linear scales. The 
open squares indicate the total current and the circles indicate the partial current of M + . The solid line is a plot of the current- 
overpotential equation (3) for io = 3.80 x 10 6 A/m 2 , iu m = —2.15 x 10 6 A/m 2 , and v — 0.777. The dashed line is a plot of 
the linear current-overpotential equation (5) for the same parameters. The dotted box indicates the points that contribute to 
Figure 2. 



• that currents in excess of the limiting current for 
the more noble species result in the deposition of 
alloys, 

• that there are changes in the double layer structure 
with current. 

As discussed in Ref. [1], the need to resolve the charge 
distribution in close proximity to the interface limits the 
size of the domain and the time spans we can model. 
Possibly, adaptive mesh techniques and implicit solution 
methods will enable us to examine larger domains and 
longer times. Nonetheless, our work here demonstrates 
that the phase field approach, using a very simple set of 
assumptions, can reproduce the rich behaviors of existing 
electrochemical theories and permit exploration of the re- 



lationship between double layer structure and interfacial 
kinetics. 
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